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1 Introduction 



O 

o 

rya ■ The most important factor for quantitative results in molecular dynamics sim- 

ulation are well developed force fields and models. In the present work, the 
C/3 ■ development of new models and the usage of force fields from the literature in 

large systems are presented. Both tasks lead to time consuming simulations 
that require massively parallel high performance computers. In the present 
work, new models for carbon dioxide [1] and cyclohexanol [2] were developed 
and a new method for the model development is introduced. Force fields and 
models for the simulation of PNIPAAm hydrogel in pure water and sodium 
fvj ' chloride solution were tested and verified [3] and used in simulations of nucle- 

ic*) \ ation processes [4, 5]. 

CNl ■ The simulations for all investigations were performed on the high performance 

computer HP XC 4000 at the Steinbuch Centre for Computing in Karlsruhe 
I/"") ■ (Germany), which is equipped with Opteron 2.6 GHz Dual Core processors. 

o 
o 

2 Development of models 

^ ■ 

^^ , Molecular models for applications in engineering are parameterized based on 

$_i ' quantum mechanical (QM) ab initio calculations and thermodynamic data. 

Both models presented here are rigid. For their usage in the engineering field, 
it is inrportant to be able to describe all thermodynamic properties of the sub- 
stance quantitatively. As the development of models is rather time consuming, 
a new and fast procedure is introduced. 
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Table 1. Parameters of the rigid carbon dioxide model, 
cc ec o"o eo Q rc-o 



nm kj-mol 1 nm kj-mol 1 



nm 



0.28137 0.10287 0.29755 0.83555 1.3589M0" 21 0.12869 



2.1 Carbon dioxide 

A molecular model for carbon dioxide was developed by optimizing the param- 
eters of the Lennard- Jones sites, the bond length and the quadrupole moment 
to experimental vapor-liquid equilibrium data. The resulting molecular model 
are listed in Table 1. After adjustment to these properties, it shows mean un- 
signed deviations to the experiment over the whole temperature range from 
triple point to critical point of 0.4 % in saturated liquid density, 1.8 % in va- 
por pressure, and 8.1 % in enthalpy of vaporization. The molecular model is 
assessed by comparing predicted thermophysical properties with experimental 
data and a reference equation of state for a large part of the fluid region. The 
average deviation for density and residual enthalpy is 4.5 and 1.7 %, respec- 
tively. The model is also capable to predict the radial distribution function, 
the second virial coefficient and transport properties, where the average devi- 
ations are 12 %. 

As an example, the description of the saturated densities by the model is plot- 
ted in Figure 1 in comparison to other models and experimental data. Here, 
the quantitatively good results can be seen. 

A more detailed discussion of the development and the results is given in 
Merker et al. [1]. 

2.2 Cyclohexanol 

For the development of the cyclohexanol model, a new procedure for param- 
eter adjustment to thermodynamic data via reduced units is introduced. The 
resulting parameters can be seen in Table 2. Compared to experimental data, 
this cyclohexanol model shows mean unsigned errors of 0.2 % in saturated 
liquid density and 3 % in vapor pressure over the whole temperature range 
from triple point to critical point. The model was used to predict the second 
virial coefficient and the transport properties, the average deviations from ex- 
perimental data are 0.1 1/mol and 25 %, respectively. 

As an example, the description of the saturated densities by the model is plot- 
ted in Figure 2 in comparison to experimental data. Here, the quantitatively 
good results can be seen. 

For further insight in the development and the new procedure for model de- 
velopment as well as the results see Merker et al. [2J. 
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Fig. 1. Saturated densities. Simulation results: o this work, □ EPM2 [6], A Vrabec 
et al. [7], □ Zhang and Duan [8, 9], A BBV [10, 11], □ Moller and Fischer [12], 
— EOS [13], ~k experimental critical point [13]. The inset is a magnified view of the 
critical point. 



3 Verification of models 

For large molecules with internal degrees of freedom, the models get more 
complex. Here, the usability of different force fields for the description of the 
swelling of Poly(N-isopropylacrylamide) (PNIPAAm) hydrogel was studied in 
water as well as in sodium chloride solutions. Different force fields for PNI- 
PAAm and models for sodium chloride from the literature were tested and 



Table 2. Coordinates and parameters of the LJ sites and the point charges in the 
principal axes system of the new molecular model for cyclohexanol. Bold characters 
indicate the represented atoms. 
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Fig. 2. Saturated densities of cyclohexanol: o model after first step of optimiza- 
tion, o model after reduced unit method, o final model, — DIPPR correlation [14], 
* experimental critical point [14]. Inset: Magnified view on the critical point. 



verified. 

The other field of investigation where models of different fluids were used, is 
the simulation of nucleation processes. Here, large systems with a large num- 
ber of molecules have to be simulated. In order simulated the nucleation rate, 
good molecular models for the used fluids are important as well. 



3.1 Swelling of poly (N-isopropylacrylamide) hydrogels 

Hydrogels are three-dimensional hydrophilic polymer networks. Their most 
characteristic property is their swelling in aqueous solutions by absorbing 
the solvent, which is influenced by various factors. Hydrogels can be used in 
many applications. Superabsorbers such as in diapers and contact lenses are 
the most common applications of hydrogels. 

The hydrogel which is studied in the present work is built up of poly(N- 
isopropylacrylamide) (PNIPAAm) cross-linked with N,N'-methylenebisacryl- 
amide (MBA). The degree of swelling in equilibrium of PNIPAAm is signifi- 
cantly influenced by many factors [3, 15, 16]. On the one hand, the swelling 
depends on the structure of the hydrogel itself, like the type of the monomer, 
but also the amount and type of cross-linker and of co-monomers. On the 
other hand, the environment conditions like temperature, type of solvent, 
solvent composition, salt concentration or pH-value of the solvent influence 
the swelling behavior. Varying these factors, the hydrogel typically shows a 
region where the hydrogel is swollen and a region where it is collapsed. In 
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between those two regions lies the region of conformation transition. The sol- 
vent composition or temperature, which is characteristic for that transition, 
is labeled here with ((9-solvent and ©-temperature). The (9-conditions are 
typically only weakly dependent on the amount of cross-linker. Therefore, the 
(9-conditions mainly depend on the environmental factors and the nature of 
the polymer chain. 

With molecular simulation it should in principle be possible to predict the 
swelling of different hydrogels upon varying any environmental factor. Two 
levels of model detail are applied for molecular simulations of hydrogels: 
coarse-grained and atomistic models. The advantage of coarse-grained mod- 
els, especially in combination with implicit solvents, is their comparatively 
low computational cost. Coarse-grained models were used to study conforma- 
tions of model hydrogels, depending on the solvent or the amount of ions in 
the solvent and the charges on the hydrogel [17, 18]. In these studies, it was 
observed that the hydrogel and the corresponding single polymer chain have 
similar conformations in good, bad and ©-solvents: Independent on whether 
a polymer chain is part of a hydrogel or a single chain, it is collapsed in bad 
solvents and stretched in good solvents and the transition occurs roughly at 
the same conditions. 

Atomistic simulations are computationally much more expensive than coarse- 
grained simulations, especially when the solvent is modeled explicitly. With 
atomistic simulations, the dynamic properties and the conformation of the 
solvent and the monomers of the hydrogels have been studied by several au- 
thors [19, 20, 21]. 

In the literature, predominantly model hydrogels were studied. Even for sim- 
ulations of real hydrogels, the validation on the basis of experimental data 
was difficult or not attempted at all. 

In this work, the swelling of PNIPAAm hydrogel was inestigated by atomistic 
molecular dynamics simulation and by experiment. The results from simula- 
tion and experiment are compared to study whether the swelling behavior of 
hydrogels can be quantitatively predicted by molecular simulation. 
From experimental data and coarse-grained simulation, it is known that 
the amount of cross-linker typically has no significant influence on the 0- 
conditions of the hydrogel [16, 17, 18, 22]. The ©-transition can be under- 
stood by studying single PNIPAAm chains as it depends on the interactions 
between the single chain and the solvent [3]. In the present work, therefore 
molecular simulations were performed for single chains in explicitly modeled 
solvent molecules for a real hydrogel. Namely, the temperature and sodium 
chloride concentration dependence of the swelling behavior of PNIPAAm in 
pure water was studied. 

For the present study, the force fields from the literature were used as pub- 
lished. No parameters were fitted. The results were compared to experimental 
data on the degree of swelling of PNIPAAm hydrogel in pure water as a func- 
tion of temperature and concentration of sodium chloride. 
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Table 3. Lennard- Jones parameters (it and e) and point charge magnitude (q e i) of 
the PNIPAAm force field OPLS [25], where e is the elementary charge. 
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Methods and molecular models 

For the molecular dynamics simulations of PNIPAAm in water, the following 
force fields from the literature were used to describe PNIPAAm: GROMOS- 
87 (G87) [23], GROMOS-96 53A6 (G53A6) [24] and OPLS-AA (OPLS) [25]. 
They were combined with the water model SPC/E [26]. For the simulations 
with salt, three different sodium chloride models were used: Chandrasekhar 
et al. [27], OPLS which uses the natrium ion model of Aqvist [28] and the 
chloride ion model of Chandrasekhar et al. [27] as well as G53A6 [24]. 
Here, only the results for the force field combination OPLS + SPC/E are 
shown, which describes the experimental data best. For results for the other 
force fields, cf. Walter et al. [3[. The Lennard- Jones and point charge param- 
eters of the PNIPAAm force field are listed in Table 3. The Lennard-Jones 
sites are characterized by the size parameter a and the energy parameter e. 
For the unlike Lennard-Jones pair interaction, a geometric mean combining 
rule for both a and e was used 

(Jij = y/CTi -aj, (1) 

eij = V e « ' e j- ( 2 ) 

For the intramolecular interactions, the method of the 1-4 interactions [29] 
was employed. 

Molecular simulations of PNIPAAm single chains were carried out with ver- 
sions 4.0.x of the GROMACS simulation package [30]. The GROMACS code 
is optimized for single processors and also for massively parallel machines. It 
was developed for the simulation of large molecules in solutions. 
Simulations at temperatures between 280 and 360 K were performed, namely 
at 290, 300, 310, 320 and 340 K. These allow to obtain the ©-temperature 
and the width of the transition region. For the study of the sodium chloride 
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influence on the swelling, simulations with concentrations of 5 and 10 mol/1 
sodium chloride in water were performed with the three models. These simu- 
lations allow to identify suitable salt models. 

Prior to these simulations, some preliminary studies were performed in which 
the different salt models were simulated in pure water. With these simulations, 
the thermophysical consistency of the models was tested. The used models for 
sodium chloride are listed in Table 4. 



Table 4. Lennard- Jones parameters (it and e) and point charge magnitude (q e i) of 
the salt models Chandrasekhar et al. [27], OPLS [27, 28] and G53A6 [24], where e 
is the elementary charge. 

model Na CI 

a J nm e / kj-mol q e i / e a j nm e / kJ-moP 1 q e i / e 

Chandrasekhar 0.189744 
OPLS 0.212645 

G53A6 0.257536 



For equilibration, the single PNIPAAm chain in water was simulated in the 
isobaric-isothermal ensemble over 1 to 5 • 10 7 time steps. The pressure was 
specified to be 1 bar and was controlled by the Berendsen barostat [31], the 
time step was 1 fs for all simulations. Newton's equations of motion were 
numerically solved with the leap-frog integrator [32] . For the long-range elec- 
trostatic interactions, particle mesh Ewald [33] with a grid spacing of 0.12 nm 
and an interpolation order of four was used. A cutoff radius of r c = 1.5 nm 
was assumed for all interactions. After equilibration, 2 to 4 • 10 7 production 
time steps were carried out with constant simulation parameters. Note that 
the production steps include the conformation transition as well as the simu- 
lation of the equilibrium. 
In order to analyze the results, the radius of gyration R g was calculated 

which characterizes the degree of stretching of the single chain, where m^ is 
the mass of site i and | |r^ 1 1 is the norm of the vector from site i to the center of 
mass of the single chain. The radius of gyration in equilibrium was calculated 
as the arithmetic mean over the last 5 • 10 6 time steps of the simulation 
together with its standard deviation. For the visualization of the results the 
open source program VMD [34] was used. 
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Temperature 

In a recent publication [3], the results for the simulations and experiments 
of the swelling of PNIPAAm hydrogel over the temperature are described in 
detail. Here, the results are only shortly summarized. 

For the successful force field combination OPLS + TIP4P, simulations at tem- 
peratures between 280 and 370 K were carried out in order to calculate the 
radius of gyration in equilibrium as a function of temperature and to use that 
information for determining the (9-temperature and the width of the transi- 
tion region. The simulations were performed starting from the same stretched 
initial configuration. The radius of gyration in equilibrium was determined 
from a time average over the last 5 ns of the run. 

Figure 3 and Table 5 present the results for the radius of gyration as a function 
of the temperature for the force field combination OPLS + SPC/E. For tem- 
peratures below 300 K, the single chain is stretched, for temperatures above 
340 K, it is collapsed. The (9-temperature is approximately 320 K. The width 
of the transition region is approximately 40 K. 

With the force field combination OPLS + SPC/E, qualitatively correct results 
were achieved. The main difference between the experimental data (cf. Figure 
4) and the prediction by molecular simulation is that the (9-temperature is 
about 15 K higher and the width of the transition region is overestimated 
by the molecular simulation. In summary, this is an unexpectedly favorable 
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Fig. 3. Radius of gyration R g of a PNIPAAm chain of 30 monomers in 14,482 water 
molecules in equilibrium over temperature T for the force field combination OPLS 
+ SPC/E. The error bars indicate the standard deviation. 
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Table 5. Radius of gyration R g of a PNIPAAm chain of 30 monomers in 14,482 
water molecules in equilibrium over temperature T for the force field combination 
OPLS + SPC/E. The numbers behind ± denote the standard deviation. 
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Fig. 4. Degree of swelling of the PNIPAAm hydrogel with different amounts of 
cross-linker as a function of temperature. Symbols: experimental data, lines: guide 
for the eye. The error bars denote the standard deviation. 



agreement between the prediction by molecular simulation and the experimen- 
tal data, especially when considering that the force fields were not trained on 
such type of data and no adjustments were made. 
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Sodium chloride concentration 

In the preliminary study, the three models for sodium chloride were simulated 
in SPC/E water at 298 K for different salt concentration. In these simula- 
tions the density of the electrolyte solution was measured. For comparison 
with experimental data [35] , the density was reduced by the density of pure 
water ad plotted over the mole fraction of sodium chloride. The results can 
be seen in Figure 5. It is obvious that the G53A6 model describes the experi- 
ment best. The results can be interpreted by looking on the water models the 
sodium chloride models were developed with: Chandrasekhar was developed 
with TIP4P [36], OPLS with TIP4P/SPC (two different models were used) 
and G53A6 with SPC [37]. 
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Fig. 5. Reduced density p/po of an aqueous sodium chloride solution as function 
of the salt concentration, given in mole fraction XNaCi for different sodium chloride 
models at 398 K in SPC/E water. The black line denotes experimental data from 
Zhang and Han [35] 



In clarify, which force field is able to predict the swelling of PNIPAAm in aque- 
ous solutions of sodium chloride, a single PNIPAAm chain of 30 monomers 
was simulated in water with 5 and 10 mol/1 of sodium chloride. The results 
as radius of gyration over the simulated time are plotted in Figure 6. The 
simulations were started throughout with a stretched initial configuration of 
the single chain. The salt models G53A6 and Chandrasekhar both yield a 
collapsed conformation in equilibrium for 10 mol/1 and a conformation in be- 
tween the stretched and the collapsed one for 5 mol/1. The model OPLS yields 
the opposite result. In experiments, the PNIPAAm hydrogel is swollen in pure 
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water and collapsed at high of sodium chloride concentrations (cf. Figure 7). 

Computational demand 

All presented simulations in Section 3.1 were carried out with the MPI based 
molecular simulation program GROMACS. The parallelization of the molec- 
ular dynamics part of GROMACS is based on the eighth shell domain decom- 
position method [30]. With GROMACS, typical simulation runs to determine 
the radius of gyration in equilibrium employ 64-128 CPUs running for 1-3 
days. For these simulations very large systems must be considered compris- 
ing typically about 58 800 interaction sites. Table 8 demonstrates the good 
scaling of the program on the HP XC4000 cluster at the Rechenzentrum of 
the Universitat Karlsruhe (TH). For these simulations a maximum memory 
of 284 MB and a maximum virtual memory of 739 MB was used. 

3.2 Nucleation of fluids 

The key properties of nucleation processes are the nucleation rate J that in- 
dicates how many clusters of the emerging phase appear in a given volume 
per time and the height Ail* of the free energy barrier that must be over- 
come to form stable nuclei. The most widespread approach for calculating 
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Fig. 6. Radius of gyration R a of a PNIPAAm chain of 30 monomers in 14,482 water 
molecules over simulation time t for the force field combination OPLS + SPC/E with 
different sodium chloride models at 300 K. 
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Fig. 7. Degree of swelling of the PNIPAAm hydrogel with different amounts of 
cross-linker as a function of sodium chloride concentration. Symbols: experimental 
data, lines: guide for the eye. 
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Fig. 8. Weak scaling of the massively parallel program GROMACS 4.0.3 on 
HP XC4000, simulation time reduced with the simulation time on one processor 
over number of processors. A PNIPAAm chain with 30 monomers in 14 482 water 
molecules was simulated. 
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these quantities is the classical nucleation theory [38], which has significant 
shortcomings, e.g., it overestimates Afl* for homogeneous vapor to liquid nu- 
cleation [39]. A more accurate theory of homogeneous nucleation, which is 
sought after, would also increase the reliability for more complex applications 
such as heterogeneous and ion-induced nucleation in the earth's atmosphere. 
An important problem is that the basic assumptions underlying the classical 
approach do not apply to nanoscopic nuclei [40] . Properties of such nuclei are 
hard to investigate experimentally, but are well accessible by molecular sim- 
ulation. For instance, equilibria [41] and vaporization processes [42] of single 
liquid droplets can be simulated to obtain the surface tension as well as heat 
and mass transfer properties of strongly curved interfaces. Similarly, very fast 
nucleation processes that occur in the immediate vicinity of the spinodal are 
experimentally inaccessible, whereas they can be studied by Monte Carlo [43] 
and molecular dynamics [44] simulation of systemswith a large number of par- 
ticles. Hence, molecular simulation is crucial for the further development of 
nucleation theory. 

In the recent investigations [4, 5], grand canonical molecular dynamics with 
McDonald's daemon is discussed and applied for sampling both nucleation 
kinetics and steady-state properties of a supersaturated vapour. Here, a se- 
ries of simulations is conducted for the truncated and shifted Lennard- Jones 
fluid which accurately describes the fluid phase coexistence of noble gases and 
methane. 



4 Conclusion 

In Section 2 it was shown, how molecular models can be developed. Two new 
models that carbon dioxide and cylcohexanol were introduced. Both models 
are able to predict many different properties they were not adjusted to quan- 
titatively correct. 

For the model development, also a new procedure via reduced units was intro- 
duced and succesfully used in the development of the model for cyclohexanol. 

In Section 3, the verification and usage of force fields and models for dif- 
ferent applications was discussed. 

The temperature dependence of swelling of PNIPAAm hydrogels in pure wa- 
ter was studied both experimentally and by molecular simulation. 
It was found that it is possible to study the conformation transition of PNI- 
PAAm hydrogel by molecular simulation of the single polymer chain. The force 
field combinations G53A6 + TIP4P and OPLS + SPC/E yield the stretched 
and the collapsed conformations as well as the conformation transition. Four 
other studied force field combinations yield only the collapsed conformation. 
With the force field combination OPLS + SPC/E, it is possible to qualita- 
tively predict the swelling of the hydrogel as a function of temperature. A 
©-temperature for the PNIPAAm single chain in pure water of approximately 
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320 K and a width of the conformation region of approximately 40 K was 
predicted. 

All three tested models for sodium chloride, namely Chandrasekhar, OPLS 
and G53A6, were able to describe the density of aqueous electrolyte solutions 
qualitatively. With the model G53A6 it was also possible to describe the den- 
sity quantitatively. Different conformations of the PNIPAAm hydrogel as a 
function of sodium chloride concentration in water could be described with 
all three models for the salt. But the OPLS force field describes a stretching 
of the chain for high sodium chloride concentration of 10 mol/1. This is the 
opposite of what the experiment shows. 

Molecular simulation is a useful tool for predicting conformation transitions 
of real hydrogels. Accurate quantitative results can be expected, if suitable 
force fields are available. 

It was also shown that it is possible to simulate the nucleation of real fluids 
and determine the nucleation rate from these simulations. Grand canonical 
molecular dynamics with McDonald's daemon was established as a method 
for steady-state simulation of nucleating vapours at high supersaturations. 
A series of simulations was conducted for the truncated and shifted Lennar- 
Jones fluid. The classical nucleation theory was found to underpredict the 
nucleation rate below the triple point. 
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